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We critically discuss and review the general ideas behind single- and multi-site coarse-grained 
(CG) models as applied to macromolecular solutions in the dilute and semi-dilute regime. We first 
consider single-site models with zero-density and density-dependent pair potentials. We highlight 
advantages and limitations of each option in reproducing the thermodynamic behavior and the 
large-scale structure of the underlying reference model. As a case study we consider solutions 
of linear homopolymers in a solvent of variable quality. Secondly, we extend the discussion to 
multi-component systems presenting, as a test case, results for mixtures of colloids and polymers. 

Specifically, we found the CG model with zero-density potentials to be unable to predict fluid- 
fluid demixing in a reasonable range of densities for mixtures of colloids and polymers of equal 
size. For larger colloids, the polymer volume fractions at which phase separation occurs are largely 
overestimated. CG models with density-dependent potentials are somewhat less accurate than 
models with zero-density potentials in reproducing the thermodynamics of the system and, although 
they presents a phase separation, they significantly underestimate the polymer volume fractions 
along the binodal. Finally, we discuss a general multi-site strategy, which is thermodynamically 
consistent and fully transferable with the number of sites, and that allows us to overcome most of 
the limitations discussed for single-site models. 

PACS numbers: 61.25.he, 65.20.De, 82.35.Lr 


I. INTRODUCTION 

Macromolecular fluids are systems characterized by a wide separation of time and length scales. Length scales range 
from the local atomic scale (« A) to the dimensions of the molecule (« 1-100 nm), to larger scales, of the order of /im, if 
mesoscopic particles are included or the system exhibits self-aggregation into supramolecular structures. Analogously, 
time scales vary from 1 ps (local atomic motion) to /is and beyond, which characterize the conformational relaxation of 
the chain. Typical examples are polymers, either synthetic and/or biological, and mixtures with other mesoparticles, 
like colloids. These soft-matter systems show complex physical behaviors, including a variety of fluid-fluid and fluid- 
solid transitions, glassy behavior, microphase separation and supramolecular self-assembly, etc. When focusing on 
these large-scale phenomena, local chemical details are often of little interest. Therefore, it is useful to develop a 
simplified description of the system without including too many details on the microscopic atomic scale. In this 
strategy, commonly referred to as coarse graining HHH, one develops models in which most of the internal degrees 
of freedom of the macromolecule are traced out, projecting the reference macromolecular system onto a system with 
only a limited number of interaction sites. For colloidal systems, one can use simple statistical-mechanics models, 
such as the simple hard-sphere model or soft-core generalizations thereof. 

A central problem in the coarse-graining approach is the derivation of the effective interactions among the reduced 
interaction sites, preserving/reproducing at the same time the thermodynamic properties and the large-scale features 
of the underlying microscopic reference system. For this purpose, a number of different strategies have been proposed in 
the last two/three decades [IH1]: structure-based methods, energy-based methods, force-matching methods, relative- 
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entropy methods, to mention the most popular. In this work we consider the structure-based route to coarse-graining 
[SHE], with the aim of determining models that accurately reproduce both the large-scale and the thermodynamic 
properties of the system under consideration in the interesting density ranges, for instance, in those in which phase 
separation occurs. Specifically, we consider single-site models, mapping each nracromolecule on a point-like particle, 
which, as we will discuss below, should be appropriate for the investigation of the thermodynamic behavior of the 
system in the low-density regime. The key advantage of such an approach is a huge decimation of degrees of freedom 
and consequently, a noticeable speed-up of numerical simulations. Moreover, single-site interaction models can be 
studied by using integral-equation methods, which represent a powerful tool for predicting local structure and thermo¬ 
dynamics of simple liquids, whose limitations and validity range are well understood m • Similar approaches have also 
been used for multi-site molecular models, such as PRISM (polymer reference interaction-site model) in the specific 
case of macromolecules hue]. They are quite successful in concentrated regimes, for instance, for polymer melts or 
polymer nanocomposites, but not very accurate in the dilute or senridilute regime, see, for instance, Refs. PH HE] for 
a discussion in the case of polymer-colloid mixtures. Another advantage of single-site CG models could be their use 
in connection with the adaptive strategies recently developed by Delle Site and coworkers |19j: see also Ref. [20] for 
the recent Hamiltonian formulation of the method. 

The first attempt to develop a simple single-site coarse-grained (CG) model for macromolecular solutions can be 
traced back to Ref. [21] , where two polymer chains with excluded volume at infinite dilution were considered. Further 
work on the issue appeared later in Refs. [22H24] . with the correct determination of the polymer-polymer effective pair 
potential. Only recently, however, has the single-site CG strategy been employed to determine the thermodynamic 
behavior of these complex systems mwm- The same strategy has also been applied to colloid-polymer solutions, 
representing polymers as monoatomic particles, interacting by means of a suitable effective potential with the colloids. 
This is a generalization of the Asakura-Oosawa model Wi¬ 
lli all CG applications, the main issue is the determination of the effective interactions among the interaction 
sites at the CG resolution, due to their inherently many-body nature. A possibility, widely used in the literature, 
is to represent the potential as a sum of pairwise contributions as derived for vanishingly small concentration of the 
macromolecules. Such a potential has a limited range of validity, being predictive only in the dilute limit, but it has 
the advantage of being properly defined, so that the standard statistical-mechanics formalism can be used without 
ambiguities to link thermodynamic properties to averages of microscopic observables in the appropriate statistical 
ensemble. 

To extend the validity of the effective interactions in a wider range of concentrations, while keeping the single-site 
model, two options are possible. One can include progressively three-body, four-body, etc., terms into the potential 
when increasing concentration. Alternatively, one can preserve the pairwise structure of the interaction, but switch to 
concentration-dependent pair potentials. The first option provides a systematic method to improve the transferability 
of the effective potential in terms of the macromolecular concentration and preserves thermodynamic consistency. 
However, it becomes rapidly unfeasible. Indeed, it requires the computation of n-body interactions by performing the 
statistical average over the internal degrees of freedom of n macromolecules while keeping their CG sites fixed in all 
possible positions. The complexity of this task obviously grows exponentially with the order n and has been attempted 
only for the lowest values of n j3UI34j . The second option is apparently easier, since it only requires the knowledge of 
suitable two-particle distribution functions as in the case of the potentials derived in the small-concentration limit. 
However, the pair correlation function in a dense system is not only the result of the direct interaction as derived 
at zero density, but it also has a contribution mediated by the other constituents present in the system. For this 
reason, the determination of the effective pair potential from the pair correlation function is not straightforward. One 
can either use numerical schemes, such as the iterative Boltzmann inversion |35j or the inverse Monte Carlo EH S3 
methods, or apply suitable methods within the integral-equation framework of liquid-state theory. For instance, in 
the case of homopolymer solutions the hypernetted-chain approximation (HNC) works extremely well and provides 
accurate pair potentials [55] . Unfortunately, in this approach the state dependence of the potentials gives rise to 
several drawbacks and inconsistencies. For instance, results depend on the ensemble one considers [391 . Moreover, 
particular care should be used in deriving the thermodynamics, as standard thermodynamic relations do no longer 
hold [551110] . 

The single-site CG strategy has been applied to linear-polymer mvnm and star-polymer HH EH EH IHHIH 
solutions in the dilute and semidilute regime, both in good solvent and in the thermal crossover toward the 6 point [31]. 
The same strategy has been also applied to study colloid-polymer solutions, generalizing the Asakura-Oosawa-Vrij 
model [501 H5] . 

In this paper we review the single-site CG strategy, comparing quantitatively results obtained by using zero-density 
and density-dependent potentials. We review the general theory, which is then applied to solutions of linear honropoly- 
mers and to mixtures of colloids and polymers in an implicit solvent. We highlight advantages and drawbacks of each 
option and discuss the delicate interplay between state-dependency and ensemble-dependency of the interactions, a 
key ingredient when discussing phase transitions as it occurs for instance in the colloid-polymer system. Finally, 
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we will mention a general multi-site strategy, which allows us to overcome the limitations of both options, therefore 
providing a truly transferrable and consistent CG model for polymer solutions in complex situations. 

The paper is organized as follows. In Sec. [TT] we discuss the general theory behind any structure-based coarse- 
graining procedure. Then, as an example, we apply it to the specific case of homopolymer solutions, discussing in 
detail the model based on a pairwise zero-density effective potential. In Sec. Ill]we introduce state-dependent pairwise 
potentials and discuss their limitations in the homopolymer case. In Sec.|W we generalize these approaches to colloid- 
polymer mixtures and present results for their phase diagram obtained by using CG models with zero-density and 
density-dependent potentials. In Sec. [V] we briefly review a recently developed multisite strategy, which provides a 
truly consistent and transferrable CG model, able to predict accurately the physics of the underlying microscopic 
system. Finally, in Sec. |VI| we draw our conclusions and perspectives. 


II. ZERO-DENSITY SINGLE-SITE COARSE-GRAINED MODELS 


A. General theory 


In this section we outline the general theory behind any coarse-graining procedure, discussing in detail the different 
approximations involved. An explicit example (homopolymers in implicit solvent) will be discussed in Sec. IIB Our 
treatment follows closely Ref. [IQ]. In the single-site CG representation, a system of N macromolecules of L units in 
a volume V is mapped onto a liquid of point-like particles, retaining only three translational degrees of freedom (for 
d = 3) per molecule. In practice, one replaces each macromolecule with a CG particle whose position R a is related 
to the positions r a ^ of the units of the macromolecule by the linear transformation 


L 

Ra AT(T’ a ^) ^ ' C]1‘ct,ii (1) 

where C; are constant coefficients that identify the representation, satisfying i c * = 1- Typical examples are the 
center-of-mass (CM) representation, in which each macromolecule is represented as a point particle located in its CM 
(Ci = 1/L, for all i), and the mid-point (MP) representation, in which the position of the effective particle coincides 
with that of the central atom (cj = &i,Ll 2 )- In the case of homopolymers, the first choice is the usual one for the linear 
topology [2U EH [271 [32], while the second choice is most common when dealing with star polymers [23 I461T48] . 

Independently of the details of the representation, in an exact mapping the effective potential (free) energy associated 
with a given CG configuration {J? Q } can be expressed as [TO] 

Veff({Ra}) = V W (N, V) + VW({R a }) + V^({R a }) + ... (2) 

where V^ n \{R a }) represents the genuine n-body contribution, defined as: 

N 

V^ n \{R a })= J2 u^(R ai ,...,R a J. (3) 

Ctl<Ot2<-"<Otn 

The zero-body term, also called volume term , is independent of the CG configuration {i? Q }. In the present case, 
in which we only reduce the number of degrees of freedom of the molecule, it only contains the free energy of a 
single isolated macromolecule, hence V) = Nv 0 with no volume dependence. In the absence of external fields, 

the one-body term is zero because of the translational invariance of the system. Except for the volume term, each 
uW(R ai ,...,R ar , ) can be expressed recursively in terms of suitable reduced distribution functions in the zero-density 
limit as, 


u^ n \R ai ,...,R a J = -iIn {R ai ,...,flaJ - (2 ’" } « (2) (R ai , R aj ) 

- u^(R ai ,R aj ,R ak )... 

- ^2 ( n-1 ’") t/" -1 ) (R a ., R aj R ai ) (4) 


where indicates the sum over all the ordered m-ples of n objects, and is the ?rth-order distribution 
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function defined as 


<7 (n) (-Rl, . . .,R n ) = (exp(-/3W int er))R 1 ,...,R„ 

/ TILi dr^e-W^{{r«A}) e -W^{{ri,i},..,{rn, i }) J]” =1 5(R* - M({r a ,J)) 

( 5 ) 

Here f7i n t ra and Wi n t er are the intramolecular and the intermolecular interaction potential, respectively. For instance, 
the two-body interaction is given by 


u^(R 1 ,R 2 ) = -^logG^(R 1 ,R 2 ), 


( 6 ) 


which shows that u^ (Rj . R 2 ) corresponds to the potential of mean force. If the system is homogeneous and isotropic, 
this potential depends only on the distance R = |i?i — Rf\ 1 hence we can simply write it as u^(R). Analogously, the 
three-body effective potential, defined in Eq. Q, is explicitly given by 


u( 3 \R u R 2 ,R 3 ) = ~ lo g g (3 \R 1 ,R 2 ,R 3 ) - v,W(R 12 ) - uW(R 13 ) - u^ 2 \R 23 ) 
1 G^iR^R^Rs) 

p og g^{Ri 2 )g^{Ri 3 )g^ 2 \R 23 y 


(7) 


with Rij = | Ri — Rj\. 

It is always possible to relate the thermodynamics in the zero-density limit to the effective interactions. Indeed, by 
using v,( 2 \r) one can compute the second virial coefficient B 2 , defined [13] ;49] by the small-density expansion of the 
pressure P, fiP = p( 1 + B 2 p +...), with p = N/V, as 


B 2 = 



( 8 ) 


Analogous relations hold for the higher-order effective potentials. For instance, the third virial coefficient B 3 can be 
expressed as [33]: 


B 3 



j dr 12 dr 13 (e-^ 3) (^-,r 13 ,r 23 ) _ ^ 
dr 12 dr 13 (e~P u ™( r ^ - l) ( e~ pu ™ 


g -/3[u (2) (ri 2 )+u (2) (ri 3 )-|-M (2) (r 2 3)] 
( r i3) _ x'j ( e -Pu (2) (r 23 ) _ _ 


(9) 


These relations show how the thermodynamic and the structural properties are intimately related and that a struc¬ 
turally consistent mapping is necessary to preserve the correct thermodynamic behavior. Therefore, any arbitrary 
truncation of the many-body series corresponds to a lack of thermodynamic consistency between the CG and the 
original model. These relations also show that virial coefficients can be used to elucidate the merits/demerits of 
any CG model in which one only considers the first few terms (typically, only the pair potential), as discussed in 
Refs. [33] [33] for a pure polymeric system and, in a more general context, in Ref. [501 . 


B. Coarse-graining homopolymer solutions 

As an example, we apply the coarse-graining methodology to a solution of homopolymers in implicit solvent. This 
is already a CG representation of a real system, since solvent degrees of freedom have been traced out and replaced 
by suitable effective monomer-monomer interactions. In general, their determination is extremely complex. However, 
a considerable simplification occurs if one only considers dilute and semidilute solutions of very long coils—this is 
the only case we consider below. In such systems the local monomer density is very small (it vanishes in the limit 
in which the length of the polymers goes to infinity) and the (osmotic) thermodynamic behavior and the large-scale 
structure (i.e., on scales that are at least of the order of the polymer size) are universal [5Tdo3| . i.e., independent of 
the microscopic details of the monomer-solvent interactions. Thus, there is no need to trace out exactly the solvent 
degrees of freedom. For instance, to study the good-solvent regime, it is enough to consider any polymer model, which 
shows local monomer-monomer repulsion. 
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Universality also significantly constrains the effective interactions. For finite values of L , the distributions Q^ 
depend on the specific polymer model under consideration, so that also the effective interactions are model dependent. 
However, in the scaling limit, i.e., for L —> oo, the adimensional combinations R g n G^ converge to distributions 

lim R 3 g n gW(R 1 ,...,R n )=gW(b 1 ,...,b n ), (10) 

L—too * 

where R g is the zero-density radius of gyration and b a = R a /R g . The distributions are universal to a large 
extent as they depend only on the quality of the solvent. For instance, the same result is obtained by considering any 
model that is appropriate to describe polymer solutions under good-solvent conditions. Because of the universality of 
the functions g( n \ also the generic n-body effective potential u < ' n \bi 1 b 2 ,..., b n ) is universal, i.e., independent of the 
polymer model, once it is expressed in terms of the adimensional vectors b a , and uniquely specified once the quality 
of the solvent has been determined. 

The first attempts to estimate the low-order terms of the above-reported many-body expansion can be traced back 
to the seminal work of Flory and Krigbaum [2T|. They showed that, in the CM representation, the effective pair 
potential ( n = 2) is approximately Gaussian with a range of the order of the radius of gyration of a single chain. 
Though the functional form of the interaction was reasonable, their mean-held treatment predicted u^ 2 \b = 0) to 
scale as L 0 2 with the length L of the polymer, hence it diverged in the scaling, infinite-length limit. Later, simple 
scaling arguments [22], renormalization-group [23: and numerical [24] calculations confirmed the overall shape of the 
interaction but found that, in the scaling limit, the potential is independent of L and it is of the order of at small 
distances. For the MP representation the potential is no longer bounded at the origin but diverges logarithmically ]!!!) 
as b — > 0 [46 :. 

A direct estimate of the effective pair potential can be obtained by determining the pair distribution function 
and using Eq. ([ 6 |. Simulation estimates of the effective pair potential, for both the CM and MP representation, 
are reported in Fig. |T] We consider Dornb-Joyce (DJ) chains [ITT] and report results for L = 2400. For the CM 
representation, an accurate parametrization of u^(b) in the scaling limit, was determined by Pelissetto and Hansen 
[32] . using a linear combination of three Gaussian functions: 

3 

W CM( & ) = J2 aieX V(~ b2 / C i)’ ( n ) 

i=1 


with ai = 0.999225, 02 = 1.1574, 03 = —0.38505, ci = 1.24051, C 2 = 0.85647, and C 3 = 0.551876. Such a parametriza¬ 
tion was obtained by fitting scaling results obtained by extrapolating athermal SAW data. This curve falls on top of 
the DJ results, see Fig. [TJ further confirming their universality. In the MP representation, the potential u'^j P (b) has 
been discussed at length in the context of star polymers. For b —> 0 it diverges logarithmically as log(l/ 6 ) jlSl (18117)TI . 
An explicit parametrization has been given by Hsu et al. ]47j (see their results for a two-arm star polymer) 


u Mp( b ) = \ lo g 


exp(— 5b 2 ) + exp(r 7 e Sb ) 


( 12 ) 


where a = 1.869, /3 = 0.815, 7 = 0.372, S = 0.405, r = 4.5. A check of these parametrizations can be obtained by 
using Eq. ([ 8 ]), which implies 


A 2 = 



(13) 


An accurate Monte Carlo estimate of the dimensionless ratio A 2 for polymers under good-solvent conditions is 5.500(3) 
[55] . If instead we use parametrizations (111 and (12 1 to compute integral (13) we obtain A 2 = 5.48 (CM) and 
A 2 = 5.51 (MP), respectively. They are quite close to the direct full-monomer (FM) estimate, confirming the accuracy 
of the two parametrizations. 

The determination of u ^ requires the computation of a triplet correlation, a function of three independent scalar 
variables, which is a cumbersome simulation task. Therefore, we only show results for polymer configurations whose 
CG sites are on the vertices of an equilateral triangle of side b. The simulations results for the DJ model are reported 
in Fig. [lj While in the MP representation the effective potential is purely attractive and diverging to — 00 as p —>■ 0 (a 
result that can be proved on general grounds SSI El]), in the CM representation it is soft, attractive at short distance 
(p < 1), and with a small repulsive tail for p > 1. For the CM representation, four-body and five-body were also 
determined pH], at least for some particular polymer configurations. In particular, at least up to n = 5, it was found 
that the generic n-body potential at small distances is positive (repulsive) for even n’s and negative (attractive) for 
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FIG. 1: Top: Effective pair potential u^ 2 \b) in the CM (left) and MP (right) representations as a function of b = r/R g . We 
report (circles) numerical DJ results and (full lines) the corresponding parametrizations (111 and ( |12| ). Bottom: Effective 
three-body potential /3u^ (b,b,b) as computed in the CM and MP representations as a function of b = r/R g . The simulation 
results (squares) are obtained by means of Monte Carlo simulations of the DJ model. The full lines are only meant to guide 
the eye. For the CM representation, DJ results agree with those obtained by using the SAW model [311 [32| . as expected on the 
basis of universality. 


odd n’s. Moreover, the strength of the n-body potential at small distances decreases for increasing n. For the MP 
representation this behavior was proved for all values of n and generic star polymers with / arms [48]: moreover, 
it was shown that the n-body potential at full overlap b\ = ... = b n = 0 decreases logarithmically, at least for star 
polymers with a large number of arms. 

The full computation of the many-body terms is very difficult and of little practical use, since numerical simulations 
of the model using three-body or higher-order interactions would be unfeasible. Indeed, the computational cost grows 
as N n , where n is the largest order included and N is the number of constituents of the system. Therefore, in most 
of the applications the many-body expansion is truncated, considering only the zero-density pair potential of mean 
force ©■ It is therefore important to quantify the accuracy of this very simplified effective model. 

As a first check of its accuracy, we consider the universal third-virial combination A 3 = B 3 /Rg, which depends 
explicitly on the three-body interaction term, see Eq. ([9]). Thus, comparison of A 3 computed in the FM model with 
the value computed in the CG model provides us with a direct estimate of the quantitative relevance of the neglected 
three-body interactions. For polymers in the scaling limit A 3 = 9.90(2) 1)5]. Using Eq. ^ with z^ 3 ) = 0 and 
parametrizations ( 11 ) and ( 12 ) we obtain instead A 3i cm = 7.85 and A 3i mp = 4.92 for the two different representations. 
The CG model underestimates the FM value A 3 = 9.90(8) in both cases: by 21% in the CM case and by 50% in the 
MP case, respectively. This shows that three-body interactions are relevant: if they are neglected, the pressure may be 
significantly underestimated. Similar conclusions are reached by directly comparing the equation of state. A simple 
estimate of Z = /3TI/p , p = N/V , for the CG system can be obtained by using the random-phase approximation 
(RPA) pS], which is expected to be accurate for systems with soft potentials and becomes exact for large densities. 
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TABLE I: Compressibility factor Z(<&) = /3Tl/p for polymers (FM) in the scaling limit [56] and for the CG model in the 
center-of-mass (CM) and in the mid-point (MP) representation. Both HNC and Monte Carlo (MC) predictions are reported. 


$ 

FM CM-MC 

CM-HNC 

MP-MC 

MP-HNC 

0.135 1.187 1.18458(1) 

1.185 

1.17869(1) 

1.182 

0.27 

1.393 1.38167(1) 

1.382 

1.36439(1) 

1.371 

0.54 

1.854 1.80067(1) 

1.803 

1.74840(1) 

1.762 

0.81 

2.371 2.23911(1) 

2.241 

2.14190(1) 

2.162 

1.09 

2.959 2.70461(1) 

2.707 

2.55534(1) 

2.582 

2.18 

5.634 4.55607(2) 

4.559 

4.18703(1) 

4.240 

4.36 

12.23 8.29709(2) 

8.303 

7.47886(3) 

7.584 


If 4> = (4:TTRgp)/3, it predicts 


Zrpa,cm = 1 + 1.71$, 

Zrpa,mp = 1 + 1-54$. (14) 

Clearly, the CG model does not capture the correct scaling of the osmotic pressure in the semidilute regime, i.e., 
z oc $ 1 - 309 [sung, underestimating the correct result. Moreover, Z depends on the chosen CG representation, a 
dependence which would be absent in the exact mapping. [118] In particular, consistently with the results for A 3 , the 
osmotic pressure in the CM representation is always larger than the MP representation estimate. 

More quantitatively, we can compare the compressibility factor Z in a wide range of densities, representative of 
both the dilute and semi-dilute regimes. In Table |T] we report FM results [56] and estimates obtained by using the 
CG models. In the latter case, we show both simulation results and estimates obtained by using integral-equation 
methods with the HNC closure (they are fully consistent in the whole density range, confirming the accuracy of the 
HNC closure for soft potentials). The CG model based on the CM representation appears to be more accurate than 
that based on the MP representation. However, both approaches show significant deviations from the correct FM 
estimates in the semidilute regime. This can be easily understood. For $ larger than 1, polymers overlap, so that 
many-body interactions, neglected in the simple CG model with pair potentials, play a relevant role. 

The results presented so far are relative to polymers under good-solvent conditions. The CG strategy can be 
extended to describe solutions in the thermal crossover region towards the 9 point, as well [44] . Indeed, universality 
holds even in this intermediate regime if properties are expressed in terms of the polymer volume fraction 4> and of the 
Zimm-Stockmayer-Fixman [621 z-variable, z <x (T — Tg)L x ^ 2 , which combines the deviation of the temperature from 
the 6 value and the chain length L , as long as logarithmic corrections (which are relevant at the 9 point) are neglected 
[5TH55] . Equivalently, and with a more direct physical meaning, the ^-variable can be replaced by the second- 
virial combination A 2 , whose functional dependence on z, A 2 (z), has been fully characterized in Ref. [ 55] . When 
approaching the good-solvent regime, we have z — > 00 and 4 . 2 ( 00 ) = 5.500(3), while, when approaching Tg , z —> 0 
and A 2 (z) ~ 47r 3 / 2 z. In Ref. [Uj it has been shown that the CG single-site model with zero-density pair potentials 
becomes accurate in an increasingly large density range when approaching the 9 point. This can be explained by 
noting the the n-th order virial coefficient scales as z n for z —> 0, so that n-body terms become increasingly less 
relevant approaching the 6 point. 

The model discussed in Ref. [H] neglects logarithmic corrections that scale as 1/lnL, hence vanish in the critical 
limit, but which may be relevant for finite values of L. Apparently, they can be neglected when considering the second 
virial coefficient or the expansion ratio that gives the variation of the radius of gyration as solvent quality changes 
[64] - t66] (see also the supplementary material of Ref. [183 f° r an extensive discussion). However, for some different 
quantities they are relevant: for instance, they determine the deviations from the ideal behavior of the equation of 
state at the 9 point in the semidilute regime and the peculiar behavior of the third virial coefficient [32, 7FFJ [58] , 
Unfortunately, it is not possible to use CG models to account also for these logarithmic corrections. Indeed, at the 9 
point the CG model with only pair potential is not thermodynamically stable [32) (69] [70]. 







III. STATE-DEPENDENT INTERACTIONS 


A. General theory 

To enlarge the density range in which CG single-site models with pairwise effective interactions can be used, one 
strategy proposed in the past is based on deriving and using state-dependent interactions. Structurally derived state- 
dependent potentials have been mostly discussed in the context of the canonical ensemble, see Refs. Eliini HD 112 and 
references therein, hence all thermodynamic and structural quantities depend on the densitv |119| p. In this approach 
the pair potential u' 2 > (r; p) is obtained by requiring the CG model to reproduce the two-point correlation function 
^( 2 )( r ; p'j a t the given value of the density. The uniqueness of such a potential is guaranteed by Henderson’s theorem 
mm- Of course, for p —>• 0 the potential u^ 2 \r\p) converges to the potential of mean force u^ 2 \r) considered before. 
The inversion of G^(r,p) to extract u^(r;p) can be performed by means of iterative procedures. For instance, one 
could use the iterative Boltzmann inversion method El- For soft potentials the HNC inversion scheme [2911551 is also 
particularly convenient. 

Although use of u^ 2 \r; p) allows one to reproduce the pair distribution function for any value of the density, 
there is no warranty that any other structural property of the underlying system- for instance, the three-body 
correlation function -is reproduced correctly. Moreover, state-dependent interactions introduce some inconsistencies 
in the calculation of standard thermodynamic properties EM EH ED- For instance, in systems with state-independent 
potentials there are two equivalent routes to the pressure. One can define it mechanically (virial pressure), as the 
force per unit area acting on the boundaries, or thermodynamically, as the derivative of the free energy with respect 
to density. In the presence of state-dependent interactions the two definitions are no longer equivalent |M EDI- 
Moreover, in the case of density-dependent potentials none of them reproduces the correct pressure of the underlying 
system, although, at least in the low-density limit, the virial expression is closer to the correct pressure than the 
thermodynamic one |.‘19l . Another problem of the approach is that effective state-dependent potentials depend on the 
ensemble in which they have been derived [323 : the equivalence of the ensembles breaks down. Therefore, different 
thermodynamic results are obtained by using CG models defined at the same state point of the underlying system 
but in different ensembles, making the computation of phase transitions and transition lines quite challenging. As a 
general message, care is needed when using state-dependent interactions to derive the thermodynamics of the original 
system and to compute free energies. In particular, one should be careful to employ state-dependent potentials only 
in the statistical ensemble in which they have been derived. 


B. Homopolymer solutions 

To elucidate the issues related to the use of state-dependent interactions, we consider again CG single-site models 
for linear homopolymers under good-solvent condition. This case, in the CM representation, has been discussed 
extensively in the past and a complete comparison between the underlying FM system and the CG model, mainly 
focused on thermodynamic, interfacial and large-scale structural properties, has been reported in Refs. [27, 29, 38J. In 
particular, Ref. |38| reports an explicit parametrization of the density-dependent pair potential obtained by matching 
the CM-CM pair distribution function for SAWs with L = 500 monomers. Since interactions are soft and the CG 
model corresponds to a monoatomic liquid, the inversion procedure was performed by using an integral-equation 
method with the HNC closure. This method requires a minimal computational effort and provides accurate estimates 
of the thermodynamic behavior. 

In Fig. [ 2 ] we report the effective pair potential u^ 2 \b; ( f>) (CM representation) for linear polymers at $ = 0,0.4,1 
obtained by using the HNC inversion procedure. The associated has been computed by FM simulations of 

the DJ model with polymers of length L = 2400. At first glance, the potentials appear to be not very sensitive to 
the polymer volume fraction. The value at full overlap increases slightly with density in the range of polymer volume 
fractions under consideration. For larger concentrations the strength of the interaction decreases again 123 [22] , as a 
consequence of the screening of the excluded-volume interaction. Moreover, the potential has a slightly longer range 
compared to the zero-density case, ensuring the correct scaling behavior of the osmotic pressure in the semi-dilute 
regime. The accuracy of the inversion can be tested by performing Monte Carlo simulations of the CG model and 
comparing the resulting pair distribution functions with those used as targets in the inversion procedure. From the 
results shown in Fig. [2] we can conclude that the HNC inversion for the CM representation is an accurate way to 
provide structurally consistent effective pair potentials. 

The results for the effective potentials reported in Ref. {35j differ somewhat from those we have determined by using 
long DJ chains, which effectively provide results in the scaling limit. The reason is that in Ref. [35[ finite-length SAWs 
(L = 500) were considered, without performing a scaling-limit extrapolation. SAW results are affected by relatively 
large scaling corrections, which increase with density (for a discussion, see Ref. [56j). even when L is of order 10 3 . In 
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FIG. 2: Left panel: Effective pair potential u ( ' 2 ' > (6; 4>) for different densities, $ = 0, $ = 0.4, <f> = 1, as obtained by HNC inversion 
(CM representation). In the inset we compare the effective potentials obtained in the scaling limit with those appropriate for 
SAWs with L = 500 sites. Right panel: Radial distribution functions between the polymer CMs as obtained by Monte Carlo 
simulations of the FM model (lines) and of the CG model with density-dependent potential (squares and circles) for 4> = 0.4 
and 4> = 1. We also report the zero-density distribution function (full line). Data for 4> = 0.4 and $ = 1 are shifted upward 
for clarity. 


the zero-density limit, scaling corrections are clearly visible in the result for the pair potential, which is somewhat 
more repulsive than the accurate expression (11), obtained by performing a proper extrapolation to the limit L —> oo. 


In particular, the value at full overlap (6 = 0) of the potential exceeds the asymptotic one by 6 %. A similar difference 
is observed for the second virial coefficient, which takes the value 6.18 if one uses the potential of Ref. )38], t° be 
compared with the value 5.50(3) obtained in the scaling limit [55]. To further test the accuracy of the potential, we 
have determined the potential at $ = 1 by using the pair distribution function obtained from simulations of the DJ 
model, finding again a discrepancy of approximately 6 % for the value at full contact, see Fig. [2] 

Now we analyze the consistency of the results obtained by using state-dependent interactions. For this purpose, we 
consider SAWs with L = 500 as our underlying system, so that we can use the effective density-dependent potentials 
reported in Ref. [3H], which apply to a large $ interval, up to $ = 2.5. Then, we determine the chemical potential 
using three different routes, that are equivalent for systems with state-independent interactions. We report results for 


Pfi = HpR 3 q ) + P^ c \ 


(15) 


which differs from the correct chemical potential by an irrelevant, model dependent constant, but which has the 
advantage of being universal. First, we consider the HNC expression for the chemical potential [751176] (HNC-route) 


/d/tHNcW = l n ( 


47T 


$ / d 3 b [h{bf - h{b)c{b) - 2 c(b)) , 


(16) 


where h(b) = g^ 2 \b; $) — 1 and the direct correlation function c(b) is related to d> and h(b) by the Ornstein-Zernike 
relation m- 


h(bi) = c(&i) + ( 3 <F/ 47 r) J db 2 c(b 2 )h(b 12 ), (17) 

where 612 = | t>i — £> 2 1 - A second possibility consists in determining first the compressibility factor Z = (3H/p by means 
of the virial expression, 




0 


db 


( 18 ) 
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TABLE II: Polymer chemical potential /3/t computed in the density-dependent CG model appropriate to describe L = 500 
SAWs (we use the parametrization of the effective pair potential reported in Ref. EH). We report results using three different 
routes (the HNC-route, the Z-route, the K-route), as discussed in the text, and the expression reported in Ref. [78] (FBD). 
The results labelled “FM scaling” are obtained by using the FM, scaling-limit equation of state reported in Ref. [56 . 



HNC-Route Z-route K-route 

FBD 

FM scaling 

0.25 

-1.98 

-2.02 

-2.03 

-1.99 

-2.11 

0.5 

-0.26 

-0.36 

-0.43 

-0.28 

-0.62 

1.0 

3.07 

2.77 

2.43 

3.06 

1.89 

1.5 

6.83 

6.14 

5.41 

6.80 

4.35 

2. 

11.00 

9.80 

8.59 

10.96 

6.90 

2.5 

15.28 

13.73 

11.98 

15.42 

9.55 


and then in computing /t as (Z-route) 


= in (^<&) + Z vir (<I>) - 1 + £ ZviT ® 1 <%• 

Finally, we consider the compressibility route (K-route), which is based on 


/3 Aic(^) = In 



/•* K(Z) - 1 

Jo t 


(19) 


( 20 ) 


with /\(<1>) given by 


K (^)- 1 = 1 + ff ( 2 ) ( 6 ; $) - 1 


( 21 ) 


For CG models with density-dependent interactions, only the K-route provides the correct chemical potential of the 
underlying model El Eg. Indeed, since the CG procedure reproduces the pair distribution function at any density, 
K( ( b), defined in Eq. (21), is the same in the CG and in the underlying model. Hence, also fifing b) defined in Eq. (20) 
is the same. Note also that, the Z-route and the K-route both require the effective potential to be computed for all 
densities smaller than the physical density of interest, hence they have a limited predictive power. 

Results are reported in Table |Tl] and in Fig. [3] It is apparent that inconsistencies between the different routes, 
well beyond the degree of inaccuracy related to the use of the HNC method, are present even in the dilute regime. 
The three different routes provide different predictions, satisfying /3/Ihnc > Pfrz > PP>k for a ll the densities under 
consideration. As a consequence, since the K-route estimate agrees with the chemical potential of the underlying 
system, the HNC-route and the Z-route both overestimate the correct chemical potential. It is interesting to observe 
that the HNC-route results are equivalent (with a small error due to HNC approximation) to those that would be 
obtained in a direct canonical Monte Carlo simulation by employing Widom insertion method m i.e., this route 
corresponds to the estimate that would be obtained in the approach referred to as passive approach in Ref. In 
other words, the HNC-route result is the one that would be obtained by using standard thermodynamic relations, 
disregarding the density dependence of the potential. As a consequence, ensemble equivalence is satisfied. If one 
performs grand-canonical simulations at chemical potential with potential u^ 2 \b] <f>), one obtains the correct 

volume fraction $. [1201 However, the fact that ensemble equivalence is satisfied is completely unrelated to the question 
whether /S/xhnc)^) is a correct estimate of the chemical potential of the underlying system. Indeed, as our results 
show, /3/iHNc(^ > ) differs significantly from the correct result. Finally, it is interesting to compare j3fiK obtained here 
(which gives the correct chemical potential for L = 500 SAWs) with the chemical potential that is obtained by using 
the equation of state of Ref. [56] . which refers to polymers in the scaling limit. The two quantities are reported in 
the inset of Fig. [3] The SAW model clearly overestimates the scaling-limit result, deviations significantly increasing 
with <E>. 
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FIG. 3: Polymer chemical potential f5ji computed in the density-dependent CG model appropriate to describe L = 500 SAWs 
(we use the parametrization of the effective pair potential reported in Ref. [38]). We report results using three different routes 
(the HNC-route, the Z-route, the K-route), as discussed in the text, and the expression reported in Ref. [75] (FBD) In the inset 
we compare the results using the K-route (which are the same as those computed directly by using L = 500 SAWs) and the 
analogous results obtained by using the FM, scaling-limit equation of state reported in Ref. [SU- 


iv. SINGLE-SITE COARSE-GRAINED MODEL FOR MULTICOMPONENT SYSTEMS 

In this section we generalize the discussion of Sec. |TT] to multicomponent systems. In particular, we focus on 
colloidal dispersions comprising large particles, colloids, usually modeled as hard spheres, and polymers in an implicit 
solvent. These systems are particularly interesting as they show a complex phase diagram which depends crucially on 
the polymer-to-colloid size ratio: for small ratios, only fluid-solid coexistence is observed, while for larger values an 
additional fluid-fluid transition is present [73HM]- Even in the absence of an explicit solvent, the computation of the 
full phase diagram is quite difficult, especially if one is interested in polymers with a large degree of polymerization. 
Therefore, CG models represent an important tool to investigate these systems. A first class of CG model is obtained 
by integrating out all polymer degrees of freedom. The resulting CG system is a one-component model of colloids 
interacting via an effective potential. Repeating the discussion of Sec. |II A[ one obtains an effective potential with 
an infinite number of many-body terms. Computationally it is unfeasible to include more that the leading, two-body 
term. However, such a truncated model is only predictive when the polymer-to-colloid size ratio is small. A less 
extreme approach consists in integrating out only the internal degrees of freedom of the polymer, representing each 
macromolecule with a monoatomic molecule, as already discussed in Sec. m After this reduction, one obtains a 
two-component system, comprising colloids and monoatomic CG polymers, which can be studied with much more 
ease than the original system. 

Two-component single-site CG models have been considered in several papers HH [S5H9I3 and also discussed in 
Ref. [IS]. Here, we shall only discuss models in which pair potentials are determined accurately by using FM data, 
in order to assess the reliability of the single-site model with pairwise interactions (other results are summarized and 
discussed in Ref. [18]). In Refs. mm we performed a careful comparison, considering both the model defined at 
zero-density and that using potentials depending on the polymer density [5S], focusing on the solvation properties of a 
single colloid in a polymer solution and on the thermodynamics in the homogeneous phase. As expected, the model is 
only accurate if q = R g /R c is less than I ( R c is the radius of the colloid). The failure of the model when polymers are 
larger than colloids can be understood physically, by noting that, when q > 1, polymers can wrap around the colloids, 
a phenomenon that cannot be modelled correctly if polymers are represented as soft spheres. Moreover, the system is 
accurate only if the polymer volume fraction = 4TrRgP P /3 (to avoid confusion with colloidal quantities, we add a 
suffix “p” to all polymer-related quantities) is less than 1, guaranteeing that the neglected three-polymer interactions 
are small. Finally, the accuracy decreases with increasing colloid volume fraction $ c = 47rR(?p c /3 (p c = N c /V is the 
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colloid density), since the relevance of the polymer-many-colloid interactions increases in this limit. 

Here we discuss the phase diagram of polymer-colloid solutions as predicted by CG single-site models. To assess 
their accuracy we need reference results to compare with. For q = 1 we will use FM results [95H95] . To the best of our 
knowledge, there are no such results for q < 1, hence we will compare our Monte Carlo estimates for the CG model 
with the binodals obtained by using the generalized free-volume theory (GFVT) P3 S3! EMSS], which is expected to 
become increasingly accurate as q decreases. 

We consider three values of q, q = 0.5, 0.8, and 1. For the CG models with zero-density and density-dependent 
interactions, we perform standard grand-canonical simulations using a recursive umbrella-sampling algorithm [lOO i 
mm . Insertions and deletions of colloids and polymers are performed by using the cluster moves introduced by Vink 
and Horbach [1021 1103] . which considerably improve the performance of the simulation. Simulation parameters are 
the fugacities z p and z c , which are normalized so that p p R 3 = z p and p c R 3 = z c for p p , p c —>• 0. Instead of z c we shall 
usually quote /3p, c = lnz c , while, as often in the literature, instead of reporting z p , we will report the volume fraction 

(r) 

of a polymer reservoir at the same value of z p . For the zero-density CG model the reservoir volume fraction can 
be obtained by inverting the corresponding equation of state ( z p = e^^ p ) which we have parametrized as: 

^ \ C 1 + 6.05117$ p + 11.60524V + 10.25884V) 1 / 2 

Z ^ p> ~ (l + 3.42865$ p ) 1 /2 

PM®?) = ln (^ $ p) +Z< ^p) ~ 1 + J o — 1 d ^- 

A. Results for q = 1 

Before studying phase separation by using the CG model, we have determined the reference binodal, using the FM 
results of Ref. [93]. Given the computational complexity of the system, the simulated chains are relatively short. 
Therefore, the results of Ref. [94] show significant corrections to scaling, which should be taken into account before 
any comparison with the CG results. The scaling-limit binodal curve can be obtained by extrapolating the data 
of Ref. [93], along the lines of the critical-point extrapolation performed in Ref. PSj. In Sec. IV.B of Ref. [IS] we 
considered the estimates of the critical points , Fc,crit(£) and < hp,crit(^)i for three systems with L = 10,33,110 and 
approximately <7 = 1, and determined the critical point in the scaling limit. We obtained [18]: 4> CiCr it(oo) ss 0.22 and 
< f > P ,crit(oo) ~ 0.62. Analogously, if $ p ln (L, 4> c ) gives the position of the binodal for the system with chains of length 
L , we fit the data to m 


( 22 ) 

(23) 


4>“ n (L, 4> c ) * 4>£ in (4> c ) + (24) 

The curve 4>p ln (4> c ) is our estimate of the scaling-limit binodal. Another possibility, although less rigorous, is to 
rescale, for each value of the length L, the finite-L binodal so as to obtain the correct critical point. In other words, 
we set 


$“ n ($ c ) = a^ in (L,6$ c ) (25) 

with 

o = b = * C ’ Cr i L \ ■ (26) 

^V,crit(-b) ^c,crit(^) 

The binodals computed with this method turn out to be essentially independent of the value of L , supporting the 
method, and quite close to that computed by direct extrapolation. The different extrapolations are reported in Fig. [4] 
together with the corresponding finite-T results. 

Once the reference binodal was determined, we considered the single-site CG model with zero-density potentials. 
We systematically increased z p and for each value of this parameter we performed several runs with different values of 
z c , covering colloid volume fractions from 0.1 to 0.35. In all cases no sign of coexistencejl22] was observed for systems 
of size V = (17.7 R g ) 3 - Of course, one might fear that systems are too small to allow us to identify a phase transition. 
Therefore, we repeated the analysis using integral-equation methods. We considered the binary system and used the 
HNC closure for all correlations. Again, no sign of phase separation was observed. 

The evidence of a wide region of stability of the homogeneous phase, well beyond the full-monomer phase boundaries, 
is surprising. Indeed, one does not expect the single-site model to be accurate if colloids and polymers have the same 
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FIG. 4: Left: binodal curves for q = 1 obtained by using the results of Ref. [331. We report the finite-L data (L = 10, L = 33, 
L = 110), the extrapolation obtained by using Eq. (241 (’’extrap”), and the binodal obtained by the simple rescaling mentioned 
in the text (’’rescal”) starting from the results with L = 110. CP is the extrapolated critical point. Right: We report 
the FM binodal (extrap), the GFVT prediction, and that obtained in Ref. [85] using polymer-density-dependent potentials 
appropriate for L = 500 SAWs (DD-SAW). We also report two points (crosses,DD-scal) belonging to the binodal obtained 
using polymer-density-dependent potentials appropriate for scaling-limit polymers. We also report the corresponding critical 
points (CP). 



FIG. 5: Partial structure factors at k = 0 as a function of for 4> c = 0.15. We report full-monomer (circles, FM) results 
obtained by using DJ chains of length L = 2400, CG results obtained by using the zero-density model [Monte Carlo results 
(triangles, CG DI-MC) and HNC results (lines, CG DI-HNC)], and by using the density-dependent model (squares, CG DD). 
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size, hence quantitative differences are not surprising. The unexpected feature is that the CG model is not even able 
to predict the qualitative behavior of the system. 

To understand why the CG model does not show phase separation, we have determined the partial structure factors 
S a fi(k) (a,/3 = p, c) and determined their limiting value for k —?■ 0. Such quantities are indeed order parameters of 
the fluid-fluid transition. We have determined these quantities for the DJ model with chains of length L = 600 and 
for the CG model for <1> C = 0.15 and several values of . For the CG model we have determined the structure factors 
both numerically, by performing Monte Carlo simulations, and by using integral-equation methods (we use the HNC 
closure) on very large systems V = ('64i? g f 3 . [l23j Results are reported in Fig. [5j For small polymer volume fractions, 
the CG and DJ results are in full agreement, but, as <F p increases, the CG model significantly underestimates the 
structure factors. At coexistence, which should occur for $ p ss 0.7-0.8, the FM estimates are S pp ( 0) ~ 4, 5 CC (0) ~ 2.5, 
which are significantly larger than the CG estimates. More precisely, for $ p = 0.76 we obtain ,S pp (0) = 4.1 and 2.3 
for the DJ and the CG model, respectively. For 5 CC (0), we obtain correspondingly 5 CC (0) = 2.2 (DJ) and 1.3 (CG). 
If we further increase <1> P , the CG results change only slightly. We obtain S pp (0) = 2.5 and 5^(0) = 1.5 for <1> P = 1.0. 
Clearly, even increasing polymer density the system appears to be unable to develop long-range correlations. 

The results for the CG model are in contrast with those of Ref. [55], which observed phase coexistence for q = 1, using 
the model with density-dependent interactions. Quantitatively, the binodal obtained in Ref. [85] differs somewhat 
from that obtained by using the FM estimates, see Fig. [4j The results for Ref. [85] refer to SAWs with L = 500 
monomers, hence one might fear that the differences between the CG and the full-monomer results are due to the 
different reference system. To clarify the issue, we have redetermined the density-dependent potential for <F P = 1 

(r) 

using scaling-limit FM data and recomputed the position of the binodal for such a value of . We find coexistence 
between (4> c ,$ p ) = (0.04,0.86) and (0.34,0.12). These two points are also reported in Fig. They show that the 
scaling-limit binodal computed by using the density-dependent potential is not very different from that computed 
by Ref. [85] and still significantly below the FM binodal. [1241 It is interesting to note that the GFVT binodal is 
essentially on top of the binodal of Ref. [85]. In view of the previous discussion, however, such an agreement looks 
accidental. 

To understand why the CG model with density-dependent potentials predicts phase separation, we have computed 
also in this case the partial structure factors. They are shown in Fig. [5] It is clear that the model with zero-density 
potential provides a better approximation to the FM results than that using density-dependent potentials. However, 
this is not relevant to obtain phase separation. The important point is that the CG model with density-dependent 
potentials overestimates significantly S pp (0) and S CC (Q), hence it exhibits phase separation, while the model with zero- 
density potentials, although more accurate in the considered range of densities, underestimates S' PP (0) and S cc ( 0), so 
that no transition occurs, at least in the range we investigated. 


B. Results for q = 0.5 and q = 0.8 


Let us now consider the behavior for q = 0.5 and 0.8. In this case we do not have reference FM results to compare 
with. Therefore, we use the GFVT predictions that are expected to become increasingly accurate as q decreases. 
Moreover, the FM results for q = 1 provide as an upper bound in <F p on the correct binodal. For a given value of <E> C , 
phase separation for q < 1 should occur at polymer volume fractions that are smaller than those at which coexistence 
occurs for q = 1. We limited our investigation here to the CG model with density independent potential. 

To identify the coexistence line, we proceed as follows. We fix z p and determine the distribution of N c and N p 
for several values of z Cl either directly or by applying the standard reweighting method 1 104] . Then, the value of z c 
corresponding to the binodal, z^ ln {z p ), is obtained by applying the usual equal-area criterion: the areas below the 
two peaks characterizing the distributions of both N c and N p should be equal. Once z^ la {z p ) has been identified, the 
averages of N c and N p over the two peaks give the number of polymers and colloids in the two phases. Results are 
reported in Tables [TTT] and |TV] for q = 0.5 and q = 0.8, respectively. They have been obtained using reasonably large 
cubic systems, of side 31.2 R g and 23.1i? g for q = 0.5 and q = 0.8, respectively. We expect size effects to be negligible, 
except possibly close to the critical point. 

To identify the critical point we use the method of Wilding m, exploiting the fact that the transition is in 
the same universality class as the three-dimensional Ising transition. In the spin system the order parameter is the 
magnetization M , whose distribution at the critical point is known quite accurately [M]: 


P(M ) = A exp 



c + a 


M 1 

K 


(27) 


with A = 0.486642, a = 0.158, c = 0.776. The normalization constant Mq can be determined by noting that 
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TABLE III: Binodal line for q = 0.5. We report the values of 4> c and <f> p in the colloid-gas (g) and in the colloid-liquid (1) 
phase. 


p 

$ (9) 

<# 

4 9) 

4° 

0.824 

0.164 

0.331 

0.336 

0.575 

0.827 

0.157 

0.340 

0.325 

0.588 

0.831 

0.149 

0.3505 

0.314 

0.603 

0.8345 

0.140 

0.361 

0.301 

0.619 

0.838 

0.132 

0.372 

0.288 

0.635 

0.842 

0.123 

0.382 

0.275 

0.651 

0.846 

0.115 

0.391 

0.2645 

0.666 

0.8495 

0.107 

0.400 

0.255 

0.680 

0.853 

0.101 

0.407 

0.246 

0.693 


TABLE IV: Binodal line for q = 0.8. We report the values of 4> c and 4>j, in the colloid-gas (g) and in the colloid-liquid (1) 
phase. 


4 ff) 4 ° 4 9) < 4 ° 


1.604 

0.177 

0.327 

0.666 

1.09 

1.610 

0.175 

0.330 

0.661 

1.10 

1.616 

0.173 

0.333 

0.655 

1.115 

1.621 

0.171 

0.337 

0.649 

1.125 

1.627 

0.168 

0.340 

0.643 

1.14 

1.633 

0.166 

0.344 

0.637 

1.15 

1.639 

0.163 

0.347 

0.631 

1.16 

1.645 

0.160 

0.351 

0.624 

1.18 

1.650 

0.157 

0.354 

0.617 

1.19 

1.656 

0.154 

0.358 

0.610 

1.20 

1.662 

0.151 

0.362 

0.604 

1.22 

1.668 

0.148 

0.365 

0.597 

1.23 

1.67346 

0.145 

0.369 

0.590 

1.24 


(M 2 ) = 0.777403,/V/g. For the mixture the order parameter analogous to the magnetization is a linear combination of 
N c and N p that can be defined as n = A(N C — aN p + b). Then, using the distributions of N c and N p computed for 
(z p ), we determine a and b by requiring (n) = 0 and the distribution to be symmetric around 


each value of z p and 2 . 


bin 


n = 0. Finally, A is determined by requiring 


= 1. Thus, for each value of the binodal we obtain a distribution 

which 


function of the variable n , which is compared with distribution (271. The best matching occurs at a value z p , 
is then identified with the critical point. The distributions at the critical point are compared with the Ising one in 
Fig. § where we report the Ising distribution with n = 1.13417M/Mo and the distributions obtained using the data 
for q = 0.5 and 0.8. The agreement is very good. The mixing of N c and N p is very small. We obtain a = 0.069, and 
0.148 for q = 0.5 and q = 0.8, respectively. This is further confirmed by the distributions of N c shown in Fig. i they 
are already quite symmetric along the binodal. 

For q = 0.5, the analysis of the data gives z PtCI it = 2.28 (equivalently = 0.823), Pn c ,c rit = 27.2. Correspond¬ 

ingly, we have 4> CjCr i t = 0.25 and <E> PiCr i t = 0.46. We have not performed a detailed analysis of the finite-box error on 
these results, but it should be of the order of 0.01 on both critical volume fractions. 


For q = 0.8 the analysis of the data gives z PjCr it = 60.11 
Correspondingly, we obtain <I> c , cr it = 0.25 and <b PiCr ; t = 0.89. 


(equivalently $^ r it = 1-621) and /?/i CjC rit = 22.9. 
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FIG. 6: Critical point distribution. The abscissa is rescaled to obtain a unit variance distribution. 



FIG. 7: Fluid-fluid binodals for q = 0.8 (left) and q = 0.5 (right). We report the single-site (CG) result and the GFVT 
prediction. For q = 0.8 we also report the binodal computed using the simplified model of Ref. [91] (CG-AP). For each binodal 
we also report the corresponding critical point (CP). 


Let us now compare the results with other estimates. For q = 0.8 it is quite evident that the single-site binodal is 
located at polymer densities that are too large. This is quite evident if we consider the location of the critical point. 
For q = 1 we estimated d> CiCr i t = 0.22 and ~ 0.62 for the full-monomer model ESI, hence the obtained estimate 

of 'Lp.crit = 0.89 is clearly far too large. It is interesting to note that the simplified model of Ref. [9Tj gives a binodal 
which is not very different from the one computed here. For q = 0.5 the single-site binodal is compatible with the 
upper bound provided by the full-monomer results with q = 1. However, comparison with the GFVT results indicate 
that, most likely, the single-site CG model predicts phase separation at values of that are too large. 








17 



FIG. 8: Distribution of <f> c at the binodal: (top) q = 0.8, (bottom) q = 0.5. We report curves for different values of &^\/3jj, c 
(they are reported in the legend). 


V. MULTI-SITE COARSE-GRAINED MODELS 

In the previous sections we have discussed single-site CG models and shown their advantages and limitations in 
describing the physics of the underlying microscopic system. Models with zero-density potentials are thermodynam¬ 
ically and structurally consistent, but are only accurate in a narrow range of densities (dilute regime) and, in the 
presence of colloids, for quite small polymer-to-colloid size ratios q. For instance, for q = 1 we have not been able 
to observe phase demixing in a reasonable range of densities, while for smaller values of q the polymer densities at 
which demixing occurs are significantly overestimated. The single-site model with state-dependent potential are only 
apparently more promising. They are not more accurate than those using zero-density potentials and moreover, as 
discussed in Sec. |III[ they are not thermodynamically consistent. 

It seems therefore quite difficult to devise an accurate and consistent CG model at the level of the single-site 
representation with pairwise interactions. On the other hand, the route of single-site models with many-body, state- 
independent interactions is also impractical. Hence, it is tempting to abandon the single-site models in favor of the 
multi-site representation. For polymeric chains in the dilute and semidilute regime, multi-site representations allow 
us to extend the density range in which they are predictive, because of the fractal nature of the chains. The volume 
occupied by a long chain in solution scales like R 3 ~ L 3l/ where v > 1/3 is the Flory scaling exponent {y = 0.5 
in 9 solvent and v ~ 0.5876 in good solvent). Moreover, since chains are self-similar objects, the volume occupied 
by each section of a chain of t monomers (blob, i 1) scales like R 3 {t) ~ £ 3 " where Rb is the radius of gyration 
of the polymer section. We have seen that the single-site CG model with pair interactions derived at zero polymer 
density provides an accurate representation of polymer solutions as far as $ p = 4ttR 3 N/3V is at most 1. Let us 
now divide each chain in n blobs of size £ in such a way that L = n£ and let us represent each polymer section by 
a single interaction site. By using state-independent interactions we can expect this model to be accurate as far as 
<I>b = AnRl(£)nN/3V = $ p n{t/L) iv = $ p n 1_3l/ < 1, i.e., for $ p < n 3u ~ l . Since 3v — 1 > 0, this relations implies that 
we can increase the density range of validity of the multi-site model and explore the semidilute regime by increasing 
the number of blobs per chain. Since interactions are derived at zero density these models preserve thermodynamic 
consistency. 

These ideas were first used in constructing simple multiblob models based on pair-wise intermolecular blob inter¬ 
actions, simple bonding intramolecular interactions and simple transferability assumptions. In Ref. m potentials 
obtained at zero-density were employed while in Ref. m the potentials were optimized to reproduce properties 
at finite density from full monomer simulations. From these first studies it was clear that multi-site models are 
more difficult than their single-site analogs. The additional difficulty comes from the fact that intramolecular and 
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FIG. 9: Grand-canonical equation of state: Polymer volume fraction <f> p versus zR 3 , where z is the fugacity defined so that 
z s=s p p for small densities. We report full-monomer (FM) results, and estimates obtained by using CG models with n = 1,4, 
and 10 sites. 


intermolecular interactions have intrinsically a many-body character. When <f> p < n 3 " -1 we can safely neglect the 
interactions among three or more polymer chains. However, the interactions among blobs belonging to the same 
molecule (intramolecular interaction) or to different molecules (intermolecular interaction) are not necessarily well 
represented as sums of pairwise potentials. Indeed, there is no obvious ’’small” physical parameter that allows us 
to expand the many-body blob potentials as a sum of two-body, three-body, ... terms. At small polymer density 
(d>p < 1), one could expect such a parameter to be the ratio between the volume of a blob and the volume of the 
chain (R^/Rg) 3 = n ~ 3v . This will indeed asymptotically control the importance of two-body nonbonding interactions 
relative to the many-body nonbonding terms. However, this parameter does not take into account the connectivity 
of the chain and its topology. It should be emphasized that the nature of the CG strategy is to replace long polymers 
by CG models with a limited number of sites. Therefore, we would like to avoid increasing the number of blobs per 
chain too much, in particular, we would like to use this parameter to control the accuracy in polymer density, not 
the accuracy of the single-chain properties. At the same time, a good CG strategy should be able to reproduce the 
structural properties of a single chain (related to the intramolecular effective potential) with any number of sites in 
a reasonable range of n. 

In Ref. [53] the first step toward such universal multi-site CG representation of polymer solutions was presented. 
We proposed to map a single linear chain in the scaling limit, the properties of which, suitably normalized, are 
universal properties of self-avoiding paths, onto a tetramer, i.e., a linear molecule with four interaction sites. The 
intramolecular force field of the tetramer was decomposed in two-body, three-body and four-body terms inspired 
by the way used in building the force fields of real molecules. Namely we introduced pair interactions between any 
pair of sites, three-body interactions in terms of bending angle potentials, and four-body interactions in terms of a 
dihedral angle potential. The potentials were obtained numerically by iteratively inverting the appropriate reduced 
probability densities. Details are given in Ref.[Si- The choice of a tetramer representation over more general n-mer 
representations with n > 4 was dictated by simplicity reasons: the tetramer is the most elaborated representation 
in which all many-body interactions (up to four body) can be parametrized as scalar functions of suitably defined 
one-dimensional variables. In a pentamer, the coupling between the two dihedral angles requires a scalar function of 
two scalar variables which is more difficult to reproduce. On the other hand, we did not limit ourself to dimer or 
trimer representations since we could not find a suitable transferability assumption allowing us to use these models 
as building blocks for more refined CG models (tetramer or more blobs) while keeping the desired accuracy. The 
tetramer CG model was found to reproduce full monomer results for polymer solutions under good solvent conditions 
up to ‘I’p ~ 2 with a ~ 5% accuracy. 

In a subsequent work [105], the tetramer model has been improved by an additional four-body contribution which 
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was found to be necessary to successfully use the tetramer model as the building block for CG representations with 
more blobs. By employing an ad hoc transferability assumption, we have shown that this new tetramer model can 
be successfully used to build CG representations with more blobs per chain, which are necessary to reach larger 
reduced polymer densities in the semi-dilute regime. It was shown that this ’’multiblob” model reproduces the leading 
correlation not explicitly considered in deriving the force held, namely the 5-body correlation between two subsequent 
dihedral angles, within ~5% of accuracy, and, more importantly, that this error does not accumulate when increasing 
the number of blobs. The result is a consistent and transferrable CG model with variable number of interaction sites 
per chain that can be used to explore the thermodynamics and the structure of polymer solutions under good solvent 
conditions. For instance, by employing a 30-blob chains, accurate results are obtained up to d> p ~ 8. In this strategy, 
all potentials are derived at zero density so that standard statistical mechanics can be applied. As an example, we 
report in Fig. [9] the equation of state obtained by grancanonical Monte Carlo simulations. Results for the single-site 
model, the tetramer model, and the decamer model are compared with the universal equation of state obtained at 
the full-monomer level in the scaling limit. Clearly, the accuracy improves systematically by increasing the number 
of sites per chain. 

The same strategy has been subsequently extended to the thermal crossover towards the 6 point JS], In this case, 
the temperature variable is expressed in terms of the parameter z (see section 2.2). As before, we have developed 
the tetramer model at selected values of z. We have found that the tetramer model is accurate in a wider range of 
reduced polymer density when approaching the 9 temperature. The transferability assumption to built multi-sites 
CG models is now more elaborated, since it must combine density and temperature, but it is equally successful, see 
Ref. [H] for details. 

More recently, we have extended the tetramer model and its n-mer extensions to colloid-polymer mixtures in a 
common good solvent [52 ]. As before we limit our intermolecular force held to pairwise central potentials. This 
strategy is successful only if many-body contributions to the free energy can be neglected. For this to happen we need 
the blob radius of gyration to be smaller than the average distance between the surface of two nearby colloids. At 

_i /o 

given reduced colloidal density d> c the average radius R s of the sphere containing a single colloid is R s /R c = <& c , 

and the average distance between the surface of two nearby colloids is d/R c = 2 (R s — R c )/R c = 2(1 — $y 3 )/$c/ 3 . 


Therefore, the condition is Rb/R c = qn u < 2(1 


$y 3 )/$y 3 or, in other terms, n > q 1 ^ 


E>!/3 


1 !/*• 


2(i-4»y 3 ). 


which 


expresses the minimum number of blobs needed for given reduce colloid density and size. Moreover, in order to treat 
a blob as a single site we need that Rb/Rc < 1 (colloidal limit) which implies n > q 1 ^. Therefore, the global condition 
on the suitable number of blobs is 


n > q 1//l/ max 


$ 


1/3 


1/v 


,2(1 -$y 3 ) 


(28) 


The term in square brackets is 1 for d? c < 0.3, while is an increasing function of <I> C for larger values. In ref. [S2| 
the homogenous phase of the mixture at various values of q was investigated. For q = 0.5 and q = 1, tetramer and 
full-monomer are found to be in full agreement up to "hp = 2. For q = 2, the tetramer slightly underestimates the 
depletion thickness. Nonetheless, it represents a significant improvement with respect to the single-blob model, which 
becomes increasingly inaccurate as <f> p increases. Also the homogeneous phase of the mixture for increasing colloid 
density was studied I92| and again the tetramer representation was found to be accurate for q = 0.5,1 at all values 
of $p, $ c in the single phase regions. At q = 2, discrepancies between the full-monomer predictions and the tetramer 
results increase with increasing both <f> p and <f> c , indicating the need of a more refined representation. By a simple 
transferability assumption for the blob-colloid potential we showed that the decamer representation can describe very 
accurately the interfacial properties of the system in the single-phase region. 


VI. CONCLUSIONS AND OUTLOOK 

In this paper we have revisited a coarse-graining strategy for polymer systems, in which a polymer chain is reduced 
to a single interaction site by tracing out the intramolecular degrees of freedom. The effective potential among the CG 
sites is inherently many-body and can be reduced to a sum of pairwise central contributions either in the low polymer- 
concentration limit or by allowing this effective pair-potential to depends on the thermodynamic state of the system. 
This CG strategy has been widely used in the past [TTl [25U2911321 [341 l47l [78l [86l l88l [89l l92l [TTn] to study polymer 
solutions in their homogeneous liquid phase and to address the theoretical study of the phase diagram of mixtures of 
non-absorbing colloids and chains of different architectures in a solvent. For homogenous polymer solutions we have 
shown the limits of validity of the CG single-site model with state-independent interactions (derived at zero polymer 
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density) and discussed the apparent improvement obtained by switching to density-dependent pair interactions. The 
latter model is indeed tuned to represent the pair correlation function at any finite polymer concentration but it 
requires the knowledge of such correlation for the underlying microscopic model, a task that need to be accomplished 
by simulating the microscopic model itself. This fact points to a limited predictive character and weakens the relevance 
of this strategy. Moreover, state-dependent interactions need to be used with care as standard thermodynamic relations 
do not hold. Equivalent routes to physical properties for state-independent potentials provide different results when 
the interaction itself depends on the thermodynamic state of the system Emm. We have explicitly discussed the 
calculation of the chain chemical potential for the homogeneous solution in Sec. |III| A main consequence of this 
inconsistency is the failure of the equivalence among different statistical ensembles even in the thermodynamic limit 


These limitation were not fully recognized at first and the state-dependent CG model was used in grand-canonical 
simulations to study the demixing transition in colloid-polymer dispersions. This CG model exhibits a demixing 


transition in qualitative agreement with experiments and phenomenological theories (GFVT). As discussed in Sec. IV 
we can now see that this agreement was accidental, since the phase transition in this model is driven by density 
fluctuations quite larger than in the underlying microscopic system. On the other hand, the CG model with zero- 
density interactions, which is thermodynamically consistent, is quite more accurate than the other model, but its 
accuracy is limited to a small range of polymer densities in the homogenous phase. Therefore, an accurate single-site 
CG model to describe the demixing transition of colloid-polymer dispersions seems to be out of reach. 

As we have briefly reviewed in the paper, in recent years ESI H09j we have developed a multi-site strategy, based on 
effective potentials derived at zero density, which is able to overcome these limitations and to provide thermodynamic 
consistent results. This new strategy is based on two main ingredients. Firstly, a suitable and transferrable repre¬ 
sentation of the intramolecular effective interaction, which allows us to keep an accurate description of single-chain 
properties for any number of sites per chain (level of coarse-graining of the model). Secondly, the possibility of using 
state-independent intermolecular interactions among CG sites of two different chains, thereby neglecting interactions 
among sites belonging to three or more chains, by adjusting the number of blobs per chain in such a way that the 
blob size be comparable to the smallest characteristic length present in the system. The strategy has been so far 
successfully applied to homopolymers under good-solvent conditions [33L 109 . and in the thermal crossover region 
towards the 6 point [2], and to colloid-polymer mixtures in the homogeneous phase E2J. We are presently using 
this approach to compute the binodal line of the colloid-polymer mixtures at the values of q for which full monomer 
data are absent. These results will provide the first quantitative determination of the phase diagram and will allow 
us to discuss on a quantitative ground the accuracy of phenomenological theories like the GFVT and the character 
of the experimental results (solvent quality). Note that it is also possible to incorporate in the model additional 
global variables, for instance the radius of gyration, as in Ref. m (see Ref. EU for the a discussion in the single-site 
context). 

The same strategy could be applied to the most disparate situations ranging from stretched and/or confined chains, 
networks, brushes, polymer nanocomposites etc. The underlying principles on which our multi-site strategy is based 
are the fractal nature of polymers and their self-similarity. On the other hand, we expect our strategy to fail if applied 
to polymer melts, since the many-body character of the effective potential cannot be truncated at the lowest order 
by just increasing the number of sites per chain. Alternative strategies have been devised for polymer melts (see the 
contribution by Guenza [112] in the same issue). However, noting that polymer-density fluctuations are much smaller 
in polymer melts than in solutions it might be possible that single-site CG models with state-dependent interaction 
provide an accurate enough description of the microscopic system and a much reduced thermodynamic inconsistency 
with respect to solutions. Other interesting systems for which only heuristic multi-site approaches have been used so 
far are di-block copolvmer solutions jll3j. grafted polvmer svstems m , and telechelic star polymer systems [ 115j . It 
would be interesting to benchmark the predictions of such heuristic model against our systematical and controllable 
strategy. 
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